
function CO2Cx32Plot(varargin)

%{
Accepts 3 optional arguments Vmax for CAII, duration of stimulus and whether to save outputs to file or not.
Units are mM and min in program, Vmax and duration specified in mM/s and mins in command line
%}

%parse inputs and assign default values if necessary

switch nargin
    case 3
        if  isnumeric(varargin{1}) && varargin{1}>=0
            Vm=varargin{1};
        else
            Vm=10;
        end
  
        if  isnumeric(varargin{2}) && varargin{2}>0 && varargin{2}<=50
            Duration=varargin{2};
        else
            Duration=10;
        end

        if varargin{3}==true || varargin{3}==false
            save=varargin{3};
        else
            save=false;
        end

    case 2
     if  isnumeric(varargin{1}) && varargin{1}>=0
            Vm=varargin{1};
     else
            Vm=10;
     end
  
     if  isnumeric(varargin{2}) && varargin{2}>0 && varargin{2}<=50
            Duration=varargin{2};
     else
            Duration=10;
     end

        save=false;

    case 1
        if  isnumeric(varargin{1}) && varargin{1}>=0
            Vm=varargin{1};
        elseif isstring(varargin{1}) && varargin{1}=="help" 
           str1="Accepts 3 optional arguments: "; 
           str1a="  Vmax for CAII (mM/s)";
           str2="   duration of stimulus (min)";
           str3="   whether to save output to csv file (true/false)";
           str4="Runs with default parameters if none specified (10,10,false)";
           disp(str1+newline+str1a+newline+str2+newline+str3+newline+str4);
           return;
        else
            Vm=10;
        end
        save=false;
        Duration=10;

    case 0
        Vm=10;
        Duration=10;
        save=false;
end


%use ode solver suitable for nonstiff equations

[t,y]=ode45(@(t,y) CO2Cx32Model(t,y,Vm,Duration),[0 100],[9.9e-3;8.5e-1;0.036;0.13;-94;0]);

subplot(7,1,1)
plot(t,y(:,1));
titleStr="Vm for CAII: " + num2str(Vm);
title(titleStr);
ylabel("Y");
%xlabel("minutes")

subplot(7,1,2)
plot(t,y(:,2));
ylabel("Z");
%xlabel("minutes")

subplot(7,1,3)
plot(t,y(:,3));
ylabel("ATP");
%xlabel("minutes");

subplot(7,1,4)
plot(t,y(:,4));
ylabel("CO2 (mM)");
%xlabel("minutes");


subplot(7,1,6)
plot(t,y(:,5));
ylabel("V (mV)");
%xlabel("minutes");

subplot(7,1,7)
plot(t,y(:,6));
ylabel("[CBF]i (mM)");
xlabel("minutes");

% Gating Cx32
HillCx=4;
Km=2.3;
CO2pow=(y(:,4)./Km).^HillCx;
gCx32=CO2pow./(1+CO2pow);


subplot(7,1,5)
plot(t,gCx32);
ylabel("gCx32");
%xlabel("minutes");


myData=[t y gCx32];

sz=size(myData);
dye1=0;

index=0;

for i=1:sz(1)
    if (myData(i,1)>=30) && (myData(i,1)<=30+Duration+2)
        dye1=myData(i,7);
        if index==0
            index=i; %this is to get the baseline value for subtraction later
        end
    end
    
end


dye1=dye1-myData(index,7);  %subtract any baseline loading to get change due to stimulation

dyeLoadTitles=["Vmax","Duration","Dye loading"];
dyeLoad=[Vm,Duration,dye1];

disp(dyeLoadTitles);
disp(dyeLoad);

if save==true
    filename="myData Vm_"+num2str(Vm)+"_"+num2str(Duration)+".csv";
    writematrix(myData,filename);
end


end



